∫Adaptive(f,a,b,N,ε) [Carlos]

bfc22725-beb8-4bf3-bbfa-80fcf15155f9 29 15 12 decimal

&[DATE] &[TIME] - &[FILENAME]

f a b N ε ∫Adaptive sum 0 : i 1 : tol i el ε : aa i el a : h i el b a - (2 / : 51 line .... fa i el a f : fc i el a h i el + f : fb i el b f : s i el h i el fa i el 4 fc i el * + fb i el + (3 / * : L i el 1 : 51 line 13 mat above => initialise the while loop i 0 > fd aa i el h i el 2 / + f : fe aa i el 3 h i el * 2 / + f : 12 mat s1 h i el fa i el 4 fd * + fc i el + 6 / (* : s2 h i el fc i el 4 fe * + fb i el + 6 / (* : 12 mat v1 aa i el : v2 fa i el : v3 fc i el : v4 fb i el : v5 h i el : v6 ε : v7 s i el : v8 L i el : 18 mat i i 1 - : s1 s2 + v7 - abs v6 < sum sum s1 s2 + (+ : v8 N ≥ Solution Insufficient Number of Intervals : i i 1 + : aa i el v1 v5 + : fa i el v3 : fc i el fe : fb i el v4 : 14 mat h i el v5 2 / : tol i el v6 2 / : s i el s2 : L i el v8 1 + : 14 mat i i 1 + : aa i el v1 : fa i el v2 : fc i el fd : fb i el v3 : h i el h i 1 - el : s i el s1 : 16 mat L i el L i 1 - el : 61 line if if 51 line while sum 41 line eval :

Example_1: singular end point

x c 1106 -^x ≤ (x 1 ≤ (& cases :

3 appVersion 0.98.6179 ≡ Timing 5 min ≡ 21 mat

x f 1 x - ln - x / :

NOTES ... 1. ∞ limit is ± ∞ at f(0) 2. 'singular' is ± ∞ at some f(?)

TOL 1013 -^:

TOL 1013 -^≤

lower limit
upper limit a TOL : b 1 TOL - : 21 sys

t0 1 time :

f a b 100 0.0000001 ∫Adaptive eval 1.644934

1 time t0 - 2 s *

t0 1 time :

S x a b 0.01 range : j 1 x rows range Y j el f a x j el 100 0.0001 ∫Adaptive eval : for x Y augment 31 line :

Zeta(2) => 1.644934066848226

x f x c * 1 10 1 1 - 2 2 mat S 1 1.645 + 40 black 1 5 mat 4 1 sys

1 time t0 - 5 s *

S S rows row 1 1.64502 12 mat

Example_2: infinite limit

infinite limit x f x 0.5 1 -^x - exp * :

'N' iterations
_ precision _iter 1000 : TOL 1015 -^: 21 sys

x f

lower limit 
upper limit a 1012 -^: b 1 : 21 sys

t0 1 time :

sanity Wolfram ... 1.49364826562485 f a b iter TOL ∫Adaptive 1.49364626562488

1 time t0 - min 1.9

π sqrt 1 erf * 1.49364826562485

ODE integration ... f(x)... discrete export

x f x ln :

end 6 :

ε 1015 -^:

x Y D x ln :

x c 10 x ≤ (x end ≤ (& cases :

sol 0 ε end 1002 x Y D dn_GearsBDF :

x f x c * sol e 2 e 1 - 2 2 mat 3.858 1.35 + 100 black 1 5 mat 4 1 sys

X sol 1 col : Y sol 2 col : 12 mat x q X Y x linterp : 21 sys

jnct x ln x q - 0 ≡ x 3.5 4 solve : 3.8579

jnct ln 1.3501

sol sol 2 col Min 2 findrows 1 1.00038 - 12 mat

a 1015 -^:

TOL 1013 -^:

f a 1100 TOL ∫Adaptive eval 1.000000000001 -

Example_3: ordinary long range ... general application

Int(f,a,b,n) ... Simpson/Romberg Kernel [Jean]

Simpson Integrator in the Romberg Transform Domain [0,1] <= [a,b] ...  
... integrates  Smath ainterp(X,Y,x), cinterp(X,Y,x), cubic(x) from  l_p_c_spline 
Integration accuracy from user command line 'n' ... NUMERIC f a b n Int w s b a - (34 / 2 w * 1 - (* 14 / 2 w * 1 - (3^* - 12 / + (* a + : w G w s f 12 w * 1 - (2^- (* : ... α 0 : β 1 : h β α - n / : 15 mat Simpson 13 / h 2 / α G β G + (* h α k h * + G k 1 n 1 - sum * + (* 23 / h * α h 2 / - k h * + G k 1 n sum * + eval : 32 / b a - (* Simpson * 41 line :

Frechet is essentially a parametric X^Y power function, not scalar wrt 'x' in Smath. Simpson-Romberg kernel resolves this inconvenience. It reduces integration on the [0 1] range then recovers. Meshing 'n' is local per application independent of WS ∫

α 1523043253 - 27354.7 - 43256.75 25388.3 - 40061 mat :

Fréchet: extreme value-II x f α 1 el α 2 el x α 4 el - (/ (α 5 el^* α 2 el x α 4 el - (/ (α 3 el^- exp * α 6 el + :

end 12 :

x c 10 x ≤ (x end ≤ (& cases :

x Y D x f :

init 0 :

f0 0 f 0 round : 403

sol init 0 end 120 x Y D dn_GearsBDF :

in-situ Simpson acc

x Frechet f 0 x 10 Int :

t0 1 time :

x f x c * sol 2 1 sys x f x c * x Frechet x c * 2 1 sys

1 time t0 - 21 s *

data X 0 end 1 range : i 1 X rows range F i el X i el Frechet : for data X F augment (: data 12 el f0 : 41 line :

data transpose 012345678910111240347018605727114941711021429243882633927642285562924529807213 mat

Gamma(x) Lanzcos Davide

z Lanczos Davde Carpi 20200514 ε 11015 -^* : p 57.156235665862923517 59.597960355475491248 - 14.136097974741747174 0.49191381609762019978 - .33994649984811888699 104 -^* .46523628927048575665 104 -^* .98374475304879564677 - (104 -^* .15808870322491248884 103 -^* .21026444172410488319 - (103 -^* .21743961811521264320 103 -^* .16431810653676389022 - (103 -^* .84418223983852743293 104 -^* .26190838401581408670 - (104 -^* .36899182659531622704 105 -^* 141 mat : g 4.7421875 : z Re 0.5 < y π π z * sin 1 z - Lanczos * / : z z 1 - : x 0.99999999999999709182 : n 1 p length range x x p n el z n + / + : for t z g + 0.5 + : y 2 π * sqrt t z 0.5 +^* t - exp * x * : 51 line if y Im abs ε ≤ y Re y if 61 line :

g=4.7421875, n=15

x Lanczos

2.5 Lanczos 1.32934038817912

2.5 Gamma 1.32934227535731

Example_4: EMERITUS Romberg-Simpson extra long range ...

x Γ x Lanczos :

x α β PDF 1 α Γ β α^* / x α 1 - (^* x β / - exp * :

x c 10 x ≤ (x 1000 ≤ (& cases :

x reduced 10 x ≤ (x 100 ≤ (& cases :

α 1.35 :

β 190 :

x α β PDF x c *

x f x α β PDF :

t0 1 time :

f 0 x 10 Int x reduced *

1 time t0 - min 0.4

Project NOTES: 1. Gamma PDF is normalized, i.e: on long range, cumulative is asymptotic to 1. 2. Decimate the quick ODE. 3. n[10] Σincremental ... plenty

t0 1 time :

The incremental cumulative companion L Δ N n Σ X 0 : Y 0 : 12 mat i 1 N range X i 1 + el L Δ i * + : X 1 el L : 21 line for X 31 line i 1 X rows range Y i el f L X i el n Int : for X Y augment 31 line :

Σincremental 01010010 Σ :

1 time t0 - min 0.1

Σincremental

x Y D x α β PDF :

init 0 :

t0 1 time :

sol_200 init 0300300 x Y D dn_GearsBDF : sol_2000 init 020002000 x Y D dn_GearsBDF : 21 sys

1 time t0 - 1.1 s *

Smath 6179 Cumulative ODE ... export

sol_200 sol_2000 2 1 sys x f x c * Gamma(x) Lanzcos Davide

z Lanczos Davde Carpi 20200514 ε 11015 -^* : p 57.156235665862923517 59.597960355475491248 - 14.136097974741747174 0.49191381609762019978 - .33994649984811888699 104 -^* .46523628927048575665 104 -^* .98374475304879564677 - (104 -^* .15808870322491248884 103 -^* .21026444172410488319 - (103 -^* .21743961811521264320 103 -^* .16431810653676389022 - (103 -^* .84418223983852743293 104 -^* .26190838401581408670 - (104 -^* .36899182659531622704 105 -^* 141 mat : g 4.7421875 : z Re 0.5 < y π π z * sin 1 z - Lanczos * / : z z 1 - : x 0.99999999999999709182 : n 1 p length range x x p n el z n + / + : for t z g + 0.5 + : y 2 π * sqrt t z 0.5 +^* t - exp * x * : 51 line if y Im abs ε ≤ y Re y if 61 line :

g=4.7421875, n=15

x Lanczos

x 2.5 i * :

2.5 Lanczos 1.32934038817912

2.5 Gamma 1.32934227535731

x Γ x Lanczos :

x α β PDF 1 α Γ β α^* / x α 1 - (^* x β / - exp * :

L H N xd U 0 : dx H L - N / : 12 mat i 1 N range U i el L dx - (dx i * + eval : for U 31 line :

ε 1015 -^:

α 0.75 :

β 250 :

x α β PDF

ε α β PDF 73

data Clear 1

L 0.1 : H 1000 : N 1000 : 13 mat

U L H N xd :

Data Σ Y 0 : vx Data 1 col : vy Data 2 col : Y 1 el 0 : i 2 vx rows range Δx i el vx i el vx i 1 - el - : Y i el Y i 1 - el 12 / vy i 1 - el vy i el + (Δx i el * (* + : 21 line for vx Y augment 61 line :

Data Y 0 : j 1 U rows range Y j el U j el α β PDF eval : for U Y augment 31 line :

Σcumul Data Σ :

U Max 999.0001

Σcumul Data Int(f,a,b,n) ... Simpson/Romberg Kernel

Simpson Integrator in the Romberg Transform Domain [0,1] <= [a,b] ...  
... integrates  Smath ainterp(X,Y,x), cinterp(X,Y,x), cubic(x) from  l_p_c_spline f a b n Int w s b a - (34 / 2 w * 1 - (* 14 / 2 w * 1 - (3^* - 12 / + (* a + : w G w s f 12 w * 1 - (2^- (* : α 0 : β 1 : h β α - n / : 14 mat Simpson 13 / h 2 / α G β G + (* h α k h * + G k 1 n 1 - sum * + (* 23 / h * α h 2 / - k h * + G k 1 n sum * + eval : 32 / b a - (* Simpson * 41 line :

x f x α β PDF :

t0 1 time :

enough 'n' on that difficult  long range n 100 :

last value cdf f L U Max n Int : 0.987

1 time t0 - 1 s *

endCumul Σcumul N row (: 999 0.989 12 mat

Yardstick integration . . .

The best know integration is the adaptive

Mathsoft proprietary. The two numerical

integrators in this document rank Smath

as equivalent to Mathsoft. How accurate

they are, can be yardsticked FD, known

and exhaustively tested very accurate wrt

to known exact integration.

]]> Σcumul U Max endCumul 1 2 el + 40 red 1 5 mat 2 1 sys

relative  ± dev'n 1 0.987 0.989 / - 0.002

x Clear 1

plug integral curve

Cumulative integration algo style ... lo mesh hi ∫ Σalgo cumullative Integrat'n algo style Y 0 : data 0 : 12 mat x lo hi mesh range : data x ∫ vectorize augment : ux data 1 col : uy data 2 col : Y 1 el 0 : i 2 ux rows range Δx i el ux i el ux i 1 - el - : Y i el Y i 1 - el 12 / uy i 1 - el uy i el + (Δx i el * (* + : 21 line for ux Y augment 91 line :

μ 0 : σ 2 : 12 mat ∫ x μ σ pdf.Normal eval : 21 sys

CDF 10 - 9.9 - 10 ∫ Σalgo :

∫ CDF 2 1 sys

x Γ x Lanczos :

α 0.75 :

β 250 :

x α β PDF 1 α Γ β α^* / x α 1 - (^* x β / - exp * :

Explain someting in there . . .

This document evaluates the cumulative integration.

The algo is the Riemann rectangle, exhaustively tested

vs more refined ancestors ... not to avail.

The algo ∫ fetches the function pdf.Normal(x,μ,σ)

CDF runs over the lo mesh, hi ,,, in fact very fast .

]]>

∫ x α β PDF eval :

lo 0.1 : mesh 1 : hi 1000 : 13 mat

t0 1 time :

CDF lo mesh hi ∫ Σalgo :

n CDF rows : 1112

CDF n row 1000 0.99 12 mat

∫ 50 ∫ * CDF 2 1 sys

1 time t0 - 8 s *

The working experimental constrained 'spline (x,v)'...m = user x v spline v 2 el x 1 m +^* 2 v 2 el * v 1 el + (x 2 m +^* - v 2 el v 1 el + (x 3 m +^* + :

v 1 - 1 - 21 mat :

m 0 :

∫ x v spline 0 x ≤ (x 1 ≤ (& * :

v 1121 mat 1 - 121 mat 11 - 21 mat 1 - 1 - 21 mat 14 mat ≡

lo 0 : mesh 0.01 : hi 1 : 13 mat

t0 1 time :

SPLINE lo mesh hi ∫ Σalgo :

n SPLINE rows : 101

known exact ... SPLINE n row 1.000 0.000 12 mat

Σalgo(lo,mesh,hi,∫) does NOT support cinterp ... Σ(Data) in lieu ...

∫ SPLINE 2 1 sys

1 time t0 - 1 s *

u b 0 u 2 - ≤ 16 / 2 u + (3^* 2 - u ≤ (u 1 - ≤ (& 16 / 2 u + (3^* 23 / 1 u + (3^* - 1 - u ≤ (u 0 ≤ (& 23 / 1 u - (3^* - 16 / 2 u - (3^* + 0 u ≤ (u 1 ≤ (& 16 / 2 u - (3^* 1 u ≤ (u 2 ≤ (& 0 cases :

∫ x b : lo 2 - : mesh 1.99 - : hi 2 : 13 mat b_spline lo mesh hi ∫ Σalgo : 31 sys

n b_spline rows : 401

b_spline n row 2112 mat

b_spline ∫ 2 1 sys

Remember this . . .

Riemann rectangular FD is almost immune to noisy

data sets. Nothing wrong to pre-smooth a bit.

]]>